/*
 * JGrass - Free Open Source Java GIS http://www.jgrass.org 
 * (C) HydroloGIS - www.hydrologis.com 
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU Lesser General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */
package eu.hydrologis.jgrass.models.h.jami;

import static eu.hydrologis.jgrass.libs.utils.JGrassConstants.isNovalue;

import java.io.File;
import java.io.PrintStream;
import java.util.ArrayList;
import java.util.Calendar;
import java.util.Collection;
import java.util.HashMap;
import java.util.Iterator;
import java.util.List;
import java.util.Map;
import java.util.Set;
import java.util.TreeMap;

import nl.alterra.openmi.sdk.backbone.ScalarSet;

import org.geotools.feature.FeatureCollection;
import org.geotools.feature.FeatureIterator;
import org.opengis.feature.simple.SimpleFeature;
import org.opengis.feature.simple.SimpleFeatureType;
import org.openmi.standard.IArgument;
import org.openmi.standard.IInputExchangeItem;
import org.openmi.standard.ILink;
import org.openmi.standard.IOutputExchangeItem;
import org.openmi.standard.ITime;
import org.openmi.standard.IValueSet;

import com.vividsolutions.jts.geom.Coordinate;
import com.vividsolutions.jts.geom.Geometry;

import eu.hydrologis.jgrass.libs.region.JGrassRegion;
import eu.hydrologis.jgrass.libs.utils.JGrassConstants;
import eu.hydrologis.jgrass.libs.utils.monitor.IProgressMonitorJGrass;
import eu.hydrologis.jgrass.libs.utils.monitor.PrintStreamProgressMonitor;
import eu.hydrologis.libs.openmi.ModelsBackbone;
import eu.hydrologis.libs.openmi.ModelsConstants;
import eu.hydrologis.libs.openmi.ModelsIOException;
import eu.hydrologis.libs.openmi.ModelsIllegalargumentException;
import eu.hydrologis.libs.utils.FluidConstants;
import eu.hydrologis.libs.utils.FluidUtils;
import eu.hydrologis.libs.utils.QuickSortAlgorithm;
import eu.hydrologis.openmi.JGrassFeatureValueSet;
import eu.hydrologis.openmi.util.HydrologisDate;

/**
 * Jami - Just another meteometeo interpolator
 * 
 * @author Silvia Franceschi (www.hydrologis.com)
 * @author Andrea Antonello (www.hydrologis.com)
 */
public class h_jami extends ModelsBackbone {

    private static final int RESOLUTION_MULTIPLIER_FOR_SIMPLIFIER = 10;

    /**
     * The type of meteo data to interpolate.
     * <p>
     * Supported types are:
     * </p>
     * <ul>
     * <li>1 - temperature</li>
     * <li>2 - pressure</li>
     * <li>3 - humidity</li>
     * <li>4 - wind</li>
     * </ul>
     */
    private int type = -1;

    /**
     * The default value for relative humidity: 70%
     */
    private double defaultRhValue = 70.0;

    /**
     * The default value for wind speed: 1 m/s
     */
    private double defaultVValue = 1.0;

    public final static String valuesID = "values";
    public final static String stationsID = "stations";
    public final static String basinsID = "basins";
    public final static String altimetryID = "altimetry";
    public final static String outputID = "out";

    private final static String modelParameters = "usage...";

    private ILink valuesLink = null;
    private ILink stationsLink = null;
    private ILink basinsLink = null;
    private ILink altimetryLink = null;
    private ILink outputLink = null;

    private IInputExchangeItem valuesInputEI = null;
    private IInputExchangeItem stationsInputEI = null;
    private IInputExchangeItem basinsInputEI = null;
    private IInputExchangeItem altimetryInputEI = null;
    private IOutputExchangeItem ouputDataOutputEI = null;

    private int stationNum = 3;
    private int stationBins = 4;
    private String stationIdField;
    private String stationElevField;
    private String basinIdField;

    private List<SimpleFeature> stationFeatures;
    private List<Coordinate> stationCoordinates;

    private List<SimpleFeature> basinFeatures;
    private List<Coordinate> basinBaricenterCoordinates;

    /**
     * The matrix holding bands elevation for every basin.
     * 
     * <p>
     * Given a number of elevation bands in every basin a matrix is created with
     * the elevation of every band. To better understand, the matrix of v values
     * is created as below (without the basin header and band column.
     * </p>
     * <table>
     * <tr>
     * <td></td>
     * <td>basin1</td>
     * <td>basin2</td>
     * <td>basin3</td>
     * <td>basin...</td>
     * <tr>
     * <td>band1</td>
     * <td>v</td>
     * <td>v</td>
     * <td>v</td>
     * <td>v...</td>
     * <tr>
     * <td>band2</td>
     * <td>v</td>
     * <td>v</td>
     * <td>v</td>
     * <td>v...</td>
     * <tr>
     * <td>band3</td>
     * <td>v</td>
     * <td>v</td>
     * <td>v</td>
     * <td>v...</td>
     * <tr>
     * <td>band...</td>
     * <td>v</td>
     * <td>v</td>
     * <td>v</td>
     * <td>v...</td>
     * </tr>
     * </table>
     */
    private double[][] bandsBasins;

    private double[] statElev;

    private double[] statId;

    /**
     * Map containing the ids of the stations and the original positions in the
     * array read from database.
     * <p>
     * This is due to the fact that the arrays are sorted at some point and we
     * need to keep track of the positions , in order to get data from unsorted
     * arrays
     * </p>
     */
    private HashMap<Integer, Integer> stationid2StationindexMap;

    /**
     * Map of basin ids and their original position in the list originally read.
     */
    private HashMap<Integer, Integer> basinid2BasinindexMap;
    private HashMap<Integer, Integer> basinindex2basinidMap;

    // private ArrayList<Geometry> basins10GeometryList;
    //
    // private ArrayList<Geometry> basins30GeometryList;
    //
    // private ArrayList<Geometry> basins80GeometryList;

    double[] radiusDistances = {10000.0, 30000.0, 80000.0};

    /**
     * Map of the list of station ids found for every elevation band of the
     * stations elevation.
     */
    private HashMap<Integer, List<Integer>> bin2StationsListMap;

    /**
     * Map of station ids and their geometry.
     */
    private HashMap<Integer, Coordinate> stationId2CoordinateMap;

    private IProgressMonitorJGrass pm;

    private JGrassRegion activeRegion;

    private int basinIdFieldIndex = -1;

    // added parameters for the DT elaboration
    private double hourForTmin = 5.0;
    private double hourForTmax = 13.0;
    private double DTd_default = 7.0;
    private double DTm_default = 7.0;
    private double tolleranceHours4Tmin = 2.0;
    private double tolleranceHours4Tmax = 2.0;

    // modificato il valore iniziale di hh_prev per le simulazioni con solo un
    // tempo...
    // rimettere il default
    // private double hh_prev = -1;
    double hh_prev = 4;

    private HydrologisDate currentTime;

    private Calendar myCalendar = Calendar.getInstance();

    private int[] cont_min_max;

    private int[] flag_Tmin;

    private int[] flag_Tmax;

    private double[] minTempPerStation;

    private double[] maxTempPerStation;

    private double[][] DTday = null;

    private double[] DTmonth = null;

    public h_jami() {
        super();
        err = FluidUtils.newPrintStream(null, System.err);
        out = FluidUtils.newPrintStream(null, System.out);
        pm = new PrintStreamProgressMonitor(out);
    }

    public h_jami( PrintStream output, PrintStream error ) {
        super();
        err = FluidUtils.newPrintStream(error, System.err);
        out = FluidUtils.newPrintStream(output, System.out);
        pm = new PrintStreamProgressMonitor(out);
    }

    public void safeInitialize( IArgument[] properties ) throws Exception {
        if (properties != null) {
            String grassDb = null;
            String location = null;
            String mapset = null;
            for( IArgument argument : properties ) {
                String key = argument.getKey();
                if (key.compareTo(ModelsConstants.GRASSDB) == 0) {
                    grassDb = argument.getValue();
                }
                if (key.compareTo(ModelsConstants.LOCATION) == 0) {
                    location = argument.getValue();
                }
                if (key.compareTo(ModelsConstants.MAPSET) == 0) {
                    mapset = argument.getValue();
                }
                if (key.compareTo("type") == 0) {
                    try {
                        type = Integer.parseInt(argument.getValue());
                    } catch (NumberFormatException e) {
                        // will handle it later
                    }
                }
                if (key.compareTo("bins") == 0) {
                    try {
                        stationBins = Integer.parseInt(argument.getValue());
                    } catch (NumberFormatException e) {
                        // will handle it later
                    }
                }
                if (key.compareTo("stationsnum") == 0) {
                    try {
                        stationNum = Integer.parseInt(argument.getValue());
                    } catch (NumberFormatException e) {
                        // will handle it later
                    }
                }
                if (key.compareTo("htmin") == 0) {
                    try {
                        hourForTmin = Integer.parseInt(argument.getValue());
                    } catch (NumberFormatException e) {
                        // will handle it later
                    }
                }
                if (key.compareTo("htmax") == 0) {
                    try {
                        hourForTmax = Integer.parseInt(argument.getValue());
                    } catch (NumberFormatException e) {
                        // will handle it later
                    }
                }
                if (key.compareTo("dtddef") == 0) {
                    try {
                        DTd_default = Integer.parseInt(argument.getValue());
                    } catch (NumberFormatException e) {
                        // will handle it later
                    }
                }
                if (key.compareTo("dtmdef") == 0) {
                    try {
                        DTm_default = Integer.parseInt(argument.getValue());
                    } catch (NumberFormatException e) {
                        // will handle it later
                    }
                }
                if (key.compareTo("tolltmin") == 0) {
                    try {
                        tolleranceHours4Tmin = Integer.parseInt(argument.getValue());
                    } catch (NumberFormatException e) {
                        // will handle it later
                    }
                }
                if (key.compareTo("tolltmax") == 0) {
                    try {
                        tolleranceHours4Tmax = Integer.parseInt(argument.getValue());
                    } catch (NumberFormatException e) {
                        // will handle it later
                    }
                }
                if (key.compareTo("stationid") == 0) {
                    stationIdField = argument.getValue();
                }
                if (key.compareTo("stationelev") == 0) {
                    stationElevField = argument.getValue();
                }
                if (key.compareTo("basinid") == 0) {
                    basinIdField = argument.getValue();
                }
            }

            if (type == -1) {
                throw new ModelsIllegalargumentException(
                        "The type of the data has not been defined. Please check your syntax.",
                        this);
            }
            if (stationIdField == null) {
                throw new ModelsIllegalargumentException(
                        "The station id field has not been supplied.", this);
            }
            if (stationElevField == null) {
                throw new ModelsIllegalargumentException(
                        "The station elevation field has not been supplied.", this);
            }
            if (basinIdField == null) {
                throw new ModelsIllegalargumentException(
                        "The basins id field has not been supplied.", this);
            }

            String activeRegionPath = grassDb + File.separator + location + File.separator + mapset
                    + File.separator + JGrassConstants.WIND;
            activeRegion = new JGrassRegion(activeRegionPath);

            valuesInputEI = ModelsConstants.createDummyInputExchangeItem(this);
            altimetryInputEI = ModelsConstants.createDummyInputExchangeItem(this);
            stationsInputEI = ModelsConstants.createFeatureCollectionInputExchangeItem(this, null);
            basinsInputEI = ModelsConstants.createFeatureCollectionInputExchangeItem(this, null);
            ouputDataOutputEI = ModelsConstants.createDummyOutputExchangeItem(this);
        }

    }

    public IValueSet safeGetValues( ITime time, String linkID ) throws Exception {
        if (outputLink == null || !linkID.equals(outputLink.getID())) {
            throw new IllegalArgumentException("A wrong output component was linked to this model.");
        }

        if (time instanceof HydrologisDate) {
            currentTime = (HydrologisDate) time;
        }
        // first ask for all you need just the first time
        if (stationCoordinates == null) {
            /*
             * get stations
             */
            out.println("Read stations data.");
            IValueSet stationsValueSet = stationsLink.getSourceComponent().getValues(time,
                    stationsLink.getID());
            FeatureCollection<SimpleFeatureType, SimpleFeature> stationFeatureCollection = ((JGrassFeatureValueSet) stationsValueSet)
                    .getFeatureCollection();
            stationCoordinates = new ArrayList<Coordinate>();
            stationFeatures = new ArrayList<SimpleFeature>();
            FeatureIterator<SimpleFeature> featureIterator = stationFeatureCollection.features();
            while( featureIterator.hasNext() ) {
                SimpleFeature feature = featureIterator.next();
                Coordinate stationCoord = ((Geometry) feature.getDefaultGeometry()).getCoordinate();
                stationCoordinates.add(stationCoord);
                stationFeatures.add(feature);
            }
            stationFeatureCollection.close(featureIterator);

            /*
             * get basins and create of every basin a buffered geometry at 10,
             * 20 and 50 km, which will be used to find the nearest stations
             * around.
             */

            IValueSet basinsValueSet = basinsLink.getSourceComponent().getValues(time,
                    basinsLink.getID());
            FeatureCollection<SimpleFeatureType, SimpleFeature> basinsFeatureCollection = ((JGrassFeatureValueSet) basinsValueSet)
                    .getFeatureCollection();
            basinBaricenterCoordinates = new ArrayList<Coordinate>();
            basinFeatures = new ArrayList<SimpleFeature>();
            featureIterator = basinsFeatureCollection.features();

            pm.beginTask("Read basins data.", basinsFeatureCollection.size());
            while( featureIterator.hasNext() ) {
                pm.worked(1);
                SimpleFeature feature = featureIterator.next();
                Coordinate baricenterCoord = ((Geometry) feature.getDefaultGeometry())
                        .getCentroid().getCoordinate();
                basinBaricenterCoordinates.add(baricenterCoord);
                basinFeatures.add(feature);
            }
            basinsFeatureCollection.close(featureIterator);
            pm.done();

            statElev = new double[stationCoordinates.size()];
            statId = new double[stationCoordinates.size()];
            stationId2CoordinateMap = new HashMap<Integer, Coordinate>();
            extractFromStationFeatures();

            stationid2StationindexMap = new HashMap<Integer, Integer>();
            for( int i = 0; i < statId.length; i++ ) {
                stationid2StationindexMap.put((int) statId[i], i);
            }
            /*
             * get the basin's id attribute index
             */
            SimpleFeature tmpfeature = basinFeatures.get(0);
            SimpleFeatureType featureType = tmpfeature.getFeatureType();
            basinIdFieldIndex = featureType.indexOf(basinIdField);
            if (basinIdFieldIndex == -1) {
                throw new IllegalArgumentException(
                        "The field of the basin id couldn't be found in the supplied basin data.");
            }
            basinid2BasinindexMap = new HashMap<Integer, Integer>();
            basinindex2basinidMap = new HashMap<Integer, Integer>();
            for( int i = 0; i < basinBaricenterCoordinates.size(); i++ ) {
                int basinid = ((Number) basinFeatures.get(i).getAttribute(basinIdFieldIndex))
                        .intValue();
                basinid2BasinindexMap.put(basinid, i);
                basinindex2basinidMap.put(i, basinid);
            }

            out.println("Creating the band's elevation for every basin matrix.");
            /*
             * get the altimetric bands matrix
             */
            IValueSet altimetryValueSet = altimetryLink.getSourceComponent().getValues(time,
                    altimetryLink.getID());
            if (altimetryValueSet instanceof ScalarSet) {
                ScalarSet altimetryScalarSet = (ScalarSet) altimetryValueSet;
                // the format is: idbasin, idfascia, elevfascia, rangefascia
                int columns = altimetryScalarSet.get(0).intValue();
                int total = altimetryScalarSet.getCount() - 1;
                int basinNum = basinBaricenterCoordinates.size();
                int fascieNum = total / basinNum / columns;

                bandsBasins = new double[fascieNum][basinNum];
                for( int i = 1; i < altimetryScalarSet.size(); i = i + 4 ) {
                    int idbasin = altimetryScalarSet.get(i).intValue();
                    int idfascia = altimetryScalarSet.get(i + 1).intValue();
                    double baricenter = altimetryScalarSet.get(i + 2);
                    Integer index = basinid2BasinindexMap.get(idbasin);
                    if (index != null)
                        bandsBasins[idfascia][index] = baricenter;
                    // TODO make it range aware
                    // double range = altimetryScalarSet.get(i + 3);
                    // bandsBasins[idfascia +
                    // 1][basinid2BasinindexMap.get(idbasin)] = baricenter
                    // + range / 2.0;
                }
            } else {
                throw new ModelsIOException(
                        "The incoming altimetry values are not scalar. Please check your syntax.",
                        this);
            }

            double[] stationBinsArrays = new double[stationBins + 1];
            double maxStatElev = statElev[statElev.length - 1];
            double minStatElev = statElev[0];
            double interval = (maxStatElev - minStatElev) / (double) stationBins;
            for( int i = 0; i < stationBinsArrays.length; i++ ) {
                stationBinsArrays[i] = minStatElev + i * interval;
            }

            /*
             * find all stations inside a elevation band for every basin
             */
            out.println("Finding all stations inside a elevation band for every basin.");
            bin2StationsListMap = new HashMap<Integer, List<Integer>>();
            for( int i = 0; i < stationBinsArrays.length - 1; i++ ) {
                List<Integer> stationsIds = new ArrayList<Integer>();
                for( int j = 0; j < statId.length; j++ ) {
                    double id = statId[j];
                    double elev = statElev[j];
                    if (elev >= stationBinsArrays[i] && elev < stationBinsArrays[i + 1]) {
                        stationsIds.add((int) id);
                    }
                }
                bin2StationsListMap.put(i, stationsIds);
            }

        }

        /*
         * get values for current timestep and order them with the stations ids
         */
        ScalarSet valuesValueSet = (ScalarSet) valuesLink.getSourceComponent().getValues(time,
                valuesLink.getID());
        double[] statValues = new double[stationCoordinates.size()];
        for( int i = 0; i < statValues.length; i++ ) {
            statValues[i] = JGrassConstants.doubleNovalue;
        }
        for( int i = 1; i < valuesValueSet.getCount(); i = i + 2 ) {
            double id = valuesValueSet.get(i);
            double value = valuesValueSet.get(i + 1);
            Integer index = stationid2StationindexMap.get((int) id);
            if (index == null)
                continue;
            statValues[index] = value;
        }

        // number of active stations for every basin
        int[] activeStationsPerBasin = new int[basinBaricenterCoordinates.size()];
        int[][] stazBacMatrix = createStationBasinsMatrix(statValues, activeStationsPerBasin);
        int[][] stations = new int[stazBacMatrix.length][stazBacMatrix[0].length];
        int contStations = 0;
        // riempimento vettori/matrici
        for( int i = 0; i < stazBacMatrix[0].length; i++ ) { // indice bacino
            contStations = 0;
            for( int j = 0; j < stazBacMatrix.length; j++ ) { // indice stazione
                if (stazBacMatrix[j][i] == 1) {
                    stations[contStations][i] = j;
                    contStations += 1;
                }
            }
        }

        ScalarSet interpolatedMeteoSet = new ScalarSet();
        int bandsNum = bandsBasins.length;
        /*
         * the number of scalars exiting at every timestep is: basin number
         * (band number + 1) where 1 is the place for the basin id
         */
        interpolatedMeteoSet.add((double) (bandsNum + 1) * basinBaricenterCoordinates.size());

        if (type == 5 || type == 6) {
            /*
             * calculate the DT month and day for each station
             */
            // System.out.println("Calculating the dayly and monthly Dt for each station...");
            rangeT(statValues);
        }

        pm.beginTask("Interpolating over bands and basins...", basinBaricenterCoordinates.size());

        // System.out.println("---");
        // for (int s = 0; s < minTempPerStation.length; s++) {
        // if (minTempPerStation[s] != 0) {
        // System.out.println("Temperature: "
        // + maxTempPerStation[s] + " "
        // + minTempPerStation[s]);
        //				
        // }
        // }

        for( int i = 0; i < basinBaricenterCoordinates.size(); i++ ) {
            pm.worked(1);
            // interpolated value for every band
            double[] interpolatedMeteoForBand = new double[bandsNum];

            int cont = 0;
            double h;
            int[] jj_av;

            // trova le stazioni che forniscono dati
            jj_av = new int[activeStationsPerBasin[i]]; // costruisco un nuovo
            // vettore jj_av con le
            // stazioni del bacino in studio
            for( int j = 0; j < activeStationsPerBasin[i]; j++ ) {
                if (type != 5 || type == 6) {
                    if (!isNovalue(statValues[stations[j][i]])) {
                        jj_av[cont] = stations[j][i]; // registro le stazioni
                        // attive
                        cont += 1;
                    }
                } else {
                    // se per la stazione j del bacino i minT e maxT sono
                    // diversi da
                    // NODATA
                    if (!isNovalue(minTempPerStation[stations[j][i]])
                            && isNovalue(maxTempPerStation[stations[j][i]])) {
                        // jj conterrà le stazioni che hanno dati di escursione
                        // termica
                        // giornaliera
                        jj_av[cont] = stations[j][i]; // registro le stazioni
                        // attive
                        cont += 1;
                    }
                }
            }

            // caso 0. se non c'e' nessuna stazione, cerco che il programma
            // sopravviva
            if (cont == 0) {
                if (type == 1) { // caso dei dati di temperatura
                    out
                            .println("ERRORE: PER IL BACINO "
                                    + i
                                    + " NON SONO DISPONIBILI DATI DI TEMPERATURA, PER QUESTO BACINO STAND-BY");
                    for( int f = 0; f < bandsNum; f++ ) { // per tutte le fasce
                        // altimetriche metto il
                        // dato a -100
                        interpolatedMeteoForBand[f] = JGrassConstants.doubleNovalue;
                    }
                } else if (type == 2) { // caso dei dati di pressione
                    out.println("  -> Per il bacino " + i
                            + " non sono disponibili dati di pressione, uso valori di default");
                    for( int f = 0; f < bandsNum; f++ ) { // per tutte le fasce
                        // altimetriche considero
                        // un'adiabatica
                        interpolatedMeteoForBand[f] = 1013.25 * Math
                                .exp(-(bandsBasins[f][i]) * 0.00013);
                    }
                } else if (type == 3) { // caso dei dati di umidità
                    out.println("  -> Per il bacino " + i
                            + " non sono disponibili dati di umidita', uso valori di default");
                    for( int f = 0; f < bandsNum; f++ ) { // per tutte le fasce
                        // altimetriche metto NODATA
                        interpolatedMeteoForBand[f] = defaultRhValue;
                    }
                } else if (type == 4) { // caso dei dati di velocità del vento
                    out
                            .println("  -> Per il bacino "
                                    + i
                                    + " non sono disponibili dati di velocita' del vento, uso valori di default");
                    for( int f = 0; f < bandsNum; f++ ) { // per tutte le fasce
                        // altimetriche metto NODATA
                        interpolatedMeteoForBand[f] = defaultVValue;
                    }
                } else if (type == 5) { // caso dei dati di escursione termica
                    // giornaliera
                    out
                            .println("  -> Per il bacino "
                                    + i
                                    + " non sono disponibili dati di escursione termica giornaliera', uso valori di default");
                    for( int f = 0; f < bandsNum; f++ ) { // per tutte le fasce
                        // altimetriche del bacino
                        // assegno all'escursione termica giornaliera il dato
                        // DTd
                        // messo nel file dei parametri
                        interpolatedMeteoForBand[f] = DTd_default;
                    }
                } else if (type == 6) { // caso dei dati di escursione termica
                    // mensile
                    out
                            .println("  -> Per il bacino "
                                    + i
                                    + " non sono disponibili dati di escursione termica mensile', uso valori di default");
                    for( int f = 0; f < bandsNum; f++ ) {
                        /*
                         *  per tutte le fasce
                         * altimetriche del bacino
                         */
                        // assegno all'escursione termica media mensile il
                        // datoDTm
                        // messo nel file dei parametri
                        interpolatedMeteoForBand[f] = DTm_default;
                    }
                }

            } else if (cont == 1) {
                // caso 1. c'e' solo una stazione presente . modello di
                // atmosfera
                // standard per T e P, valori costanti per RH e V
                for( int f = 0; f < bandsNum; f++ ) { // ciclo sulle fascie
                    // altimetriche
                    if (type == 1) { // trasformo la temp in K e calcolo T con
                        // l'adiabatica semplice
                        interpolatedMeteoForBand[f] = (statValues[jj_av[0]] + FluidConstants.tk)
                                * Math.exp(-(bandsBasins[f][i] - statElev[jj_av[0]])
                                        * FluidConstants.GAMMA
                                        / (statValues[jj_av[0]] + FluidConstants.tk))
                                - FluidConstants.tk;
                    } else if (type == 2) { // calcolo P con il gradiente
                        // adiabatico
                        interpolatedMeteoForBand[f] = statValues[jj_av[0]]
                                * Math.exp(-(bandsBasins[f][i] - statElev[jj_av[0]]) * 0.00013);
                    } else if (type == 5) {
                        // se ho una sola stazione assegno il valore della
                        // stazione a tutto il
                        // bacino
                        // altimetriche del bacino assegno il valore di
                        // escursione massima
                        // giornaliera
                        interpolatedMeteoForBand[f] = maxTempPerStation[jj_av[0]]
                                - minTempPerStation[jj_av[0]];
                        if ((maxTempPerStation[jj_av[0]] - minTempPerStation[jj_av[0]]) <= 0) {
                            interpolatedMeteoForBand[f] = DTd_default;
                        }
                    } else if (type == 6) {
                        // se ho una sola stazione assegno il valore della
                        // stazione a tutto il
                        // bacino
                        // altimetriche del bacino assegno il valore di
                        // escursione massima mensile
                        interpolatedMeteoForBand[f] = DTmonth[jj_av[0]];
                    } else { // RH e V sono costanti al variare delle fasce
                        // altimetriche
                        interpolatedMeteoForBand[f] = statValues[jj_av[0]];
                    }
                }
            } else {
                // caso 2. ci sono almeno 2 stazioni (a quote inferiori alla
                // stazioni piu' bassa considero atmosfera standard come a quote
                // superiori alla staz. piu' alta, in mezzo calcolo LAPSE RATE)
                // alloca L (vettore di dimensioni numero di stazioni attive-1)
                double[] lapseRate = new double[cont - 1];

                for( int j = 0; j < cont - 1; j++ ) { // le stazioni sono in
                    // ordine di
                    // quota
                    // L[j] e' il lapse rate tra la stazione j e j+1, puo'
                    // essere
                    // calcolato dai dati per j che va da 1 a n-1, dove n e' il
                    // numero di stazioni (cont)
                    lapseRate[j] = (statValues[jj_av[j]] - statValues[jj_av[j + 1]])
                            / (statElev[jj_av[j + 1]] - statElev[jj_av[j]]);
                }

                for( int f = 0; f < bandsNum; f++ ) { // ciclo sulle fascie
                    // altimetriche

                    // per le fasce altimetriche con quote piu' basse della
                    // quota
                    // della stazione piu' bassa prendo i dati della stazione
                    // più bassa
                    if (bandsBasins[f][i] <= statElev[jj_av[0]]) {
                        if (type == 1) { // T
                            interpolatedMeteoForBand[f] = statValues[jj_av[0]]
                                    - FluidConstants.GAMMA
                                    * (bandsBasins[f][i] - statElev[jj_av[0]]);
                        } else if (type == 2) { // P
                            interpolatedMeteoForBand[f] = statValues[jj_av[0]]
                                    - (statValues[jj_av[0]] * 0.00013)
                                    * (bandsBasins[f][i] - statElev[jj_av[0]]);
                        } else if (type == 5) {
                            interpolatedMeteoForBand[f] = maxTempPerStation[jj_av[0]]
                                    - minTempPerStation[jj_av[0]];
                            if ((maxTempPerStation[jj_av[0]] - minTempPerStation[jj_av[0]]) <= 0) {
                                interpolatedMeteoForBand[f] = DTd_default;
                            }
                        } else if (type == 6) {
                            interpolatedMeteoForBand[f] = DTmonth[jj_av[0]];
                        } else { // RH e V
                            interpolatedMeteoForBand[f] = statValues[jj_av[0]];
                        }

                        // per le fasce altimetriche con quote piu' alte della
                        // quota
                        // della stazione piu' alta prendo i dati della stazione
                        // più alta
                    } else if (bandsBasins[f][i] >= statElev[jj_av[cont - 1]]) {
                        if (type == 1) { // T
                            interpolatedMeteoForBand[f] = statValues[jj_av[cont - 1]]
                                    - FluidConstants.GAMMA
                                    * (bandsBasins[f][i] - statElev[jj_av[cont - 1]]);
                        } else if (type == 2) { // P
                            interpolatedMeteoForBand[f] = statValues[jj_av[cont - 1]]
                                    - (statValues[jj_av[cont - 1]] * 0.00013)
                                    * (bandsBasins[f][i] - statElev[jj_av[cont - 1]]);
                        } else if (type == 5) {
                            interpolatedMeteoForBand[f] = maxTempPerStation[jj_av[cont - 1]]
                                    - minTempPerStation[jj_av[cont - 1]];
                            if ((maxTempPerStation[jj_av[0]] - minTempPerStation[jj_av[0]]) <= 0) {
                                interpolatedMeteoForBand[f] = DTd_default;
                            }
                        } else if (type == 6) {
                            interpolatedMeteoForBand[f] = DTmonth[jj_av[cont - 1]];
                        } else { // RH e V
                            interpolatedMeteoForBand[f] = statValues[jj_av[cont - 1]];
                        }

                    } else {
                        int k = cont - 1;
                        if (type == 5) {
                            // per le fasce altimetriche intermedie devo
                            // interpolare tra la min e
                            // la max delle stazioni
                            do {
                                k -= 1;
                                h = statElev[jj_av[k]];
                            } while( bandsBasins[f][i] <= h );

                            // for (int j = 0; j < cont; j++) {
                            // if (f ==0 && i == 100) {
                            // System.out.println(j + " "+ statElev[jj_av[j]]);
                            // }
                            // }

                            // interpolatedMeteoForBand[f] =
                            // ((maxTempPerStation[jj_av[k]] -
                            // minTempPerStation[jj_av[k]])
                            // * (statElev[jj_av[k + 1]] - bandsBasins[f][i]) +
                            // (maxTempPerStation[jj_av[k + 1]] -
                            // minTempPerStation[jj_av[k + 1]])
                            // * (bandsBasins[f][i] - statElev[jj_av[k]]))
                            // / (statElev[jj_av[k + 1]] - statElev[jj_av[k]]);
                            interpolatedMeteoForBand[f] = ((maxTempPerStation[jj_av[k + 1]] - minTempPerStation[jj_av[k + 1]]) - (maxTempPerStation[jj_av[k]] - minTempPerStation[jj_av[k]]))
                                    * (bandsBasins[f][i] - statElev[jj_av[k]])
                                    / (statElev[jj_av[k + 1]] - statElev[jj_av[k]])
                                    + (maxTempPerStation[jj_av[k]] - minTempPerStation[jj_av[k]]);
                            // if (i == 100) {
                            // System.out.println("Banda " + f + " "
                            // + bandsBasins[f][i]);
                            // System.out.println("stazione1 " + k);
                            // System.out.println("elevazione: "
                            // + statElev[jj_av[k]]);
                            // System.out.println("stazione2 " + k + 1);
                            // System.out.println("max: "
                            // + maxTempPerStation[jj_av[k + 1]]);
                            // System.out.println("min: "
                            // + minTempPerStation[jj_av[k + 1]]);
                            // // System.out.println(statElev[jj_av[k + 1]]);
                            // }

                            if (interpolatedMeteoForBand[f] <= 0) {
                                interpolatedMeteoForBand[f] = DTd_default;
                            }
                        } else if (type == 6) {
                            // per le fasce altimetriche intermedie devo
                            // interpolare tra la min e
                            // la max delle stazioni
                            do {
                                k -= 1;
                                h = statElev[jj_av[k]];
                            } while( bandsBasins[f][i] <= h );
                            interpolatedMeteoForBand[f] = (DTmonth[jj_av[k]]
                                    * (statElev[jj_av[k + 1]] - bandsBasins[f][i]) + DTmonth[jj_av[k + 1]]
                                    * (bandsBasins[f][i] - statElev[jj_av[k]]))
                                    / (statElev[jj_av[k + 1]] - statElev[jj_av[k]]);
                        } else {
                            do {
                                k -= 1;
                                h = statElev[jj_av[k]];
                            } while( bandsBasins[f][i] <= h );
                            interpolatedMeteoForBand[f] = statValues[jj_av[k]] - lapseRate[k]
                                    * (bandsBasins[f][i] - statElev[jj_av[k]]);
                        }
                    }

                }

                // ADDED
                // controllo su RH>100 e v=0
                if (type == 3) { // RH
                    for( int f = 0; f < bandsNum; f++ ) {
                        if (interpolatedMeteoForBand[f] > 100)
                            interpolatedMeteoForBand[f] = 100;
                        if (interpolatedMeteoForBand[f] < 5)
                            interpolatedMeteoForBand[f] = 5;
                    }
                } else if (type == 4) { // V
                    for( int f = 0; f < bandsNum; f++ ) {
                        if (interpolatedMeteoForBand[f] < 0.01)
                            interpolatedMeteoForBand[f] = 0.01;
                    }
                }

                int basinid = ((Number) basinFeatures.get(i).getAttribute(basinIdFieldIndex))
                        .intValue();
                interpolatedMeteoSet.add((double) basinid);
                for( double value : interpolatedMeteoForBand ) {
                    interpolatedMeteoSet.add(value);
                }
            }
        }
        pm.done();

        return interpolatedMeteoSet;
    }

    /**
     * Creates the stations per basins matrix.
     * 
     * <p>
     * The matrix is a bitmap of the stations that will be used for every basin,
     * following the schema:
     * </p>
     * <table>
     * <tr>
     * <td></td>
     * <td>basin1</td>
     * <td>basin2</td>
     * <td>basin3</td>
     * <td>basin...</td>
     * <tr>
     * <td>station1</td>
     * <td>1</td>
     * <td>0</td>
     * <td>1</td>
     * <td>0...</td>
     * <tr>
     * <td>station2</td>
     * <td>1</td>
     * <td>0</td>
     * <td>1</td>
     * <td>0...</td>
     * <tr>
     * <td>station3</td>
     * <td>0</td>
     * <td>1</td>
     * <td>1</td>
     * <td>0...</td>
     * <tr>
     * <td>station...</td>
     * <td>0</td>
     * <td>0</td>
     * <td>0</td>
     * <td>0...</td>
     * </tr>
     * </table>
     * <p>
     * In the above case basin1 will use station1 and station2, while basin2
     * will use only station3, and so on.
     * </p>
     * 
     * @param statValues
     *            the station data values, properly ordered, containing
     *            {@link JGrassConstants#defaultNovalue novalues}.
     * @param activeStationsPerBasin
     *            an array to be filled with the number of active stations per
     *            basin.
     * @return the matrix of active stations for every basin.
     */
    private int[][] createStationBasinsMatrix( double[] statValues, int[] activeStationsPerBasin ) {
        int[][] stationsBasins = new int[stationCoordinates.size()][basinBaricenterCoordinates
                .size()];
        Set<Integer> bandsIdSet = bin2StationsListMap.keySet();
        Integer[] bandsIdArray = (Integer[]) bandsIdSet.toArray(new Integer[bandsIdSet.size()]);
        // for every basin
        for( int i = 0; i < basinBaricenterCoordinates.size(); i++ ) {
            Coordinate basinBaricenterCoordinate = basinBaricenterCoordinates.get(i);
            // for every stations band
            int activeStationsForThisBasin = 0;
            for( int j = 0; j < bandsIdArray.length; j++ ) {
                int bandId = bandsIdArray[j];
                List<Integer> stationIdsForBand = bin2StationsListMap.get(bandId);

                /*
                 * search for the nearest stations that have values.
                 */
                List<Integer> stationsToUse = extractStationsToUse(basinBaricenterCoordinate,
                        stationIdsForBand, stationId2CoordinateMap, statValues,
                        stationid2StationindexMap);
                if (stationsToUse.size() < stationNum) {
                    err.println("Found only " + stationsToUse.size() + " for basin "
                            + basinindex2basinidMap.get(i) + " and bandid " + bandId + ".");
                }

                /*
                 * now we have the list of stations to use. With this list we
                 * need to enable (1) the proper matrix positions inside the
                 * stations-basins matrix.
                 */
                // i is the column (basin) index
                // the station id index can be taken from the idStationIndexMap
                for( Integer stationIdToEnable : stationsToUse ) {
                    int stIndex = stationid2StationindexMap.get(stationIdToEnable);
                    stationsBasins[stIndex][i] = 1;
                }
                activeStationsForThisBasin = activeStationsForThisBasin + stationsToUse.size();
            }

            activeStationsPerBasin[i] = activeStationsForThisBasin;
        }

        return stationsBasins;
    }
    /**
     * @param basinBaricenterCoordinate the basin baricenter coordinate used to 
     *                      Calculate distances with stations.
     * @param stationIdsToSearch the list of stations.
     * @param stationId2CoordinateMap
     * @param statValues the array of data to consider. Used to identify stations
     *                      that have no data for an instant.
     * @param idStationIndexMap
     * @return the list of needed nearest stations.
     */
    private List<Integer> extractStationsToUse( Coordinate basinBaricenterCoordinate,
            List<Integer> stationIdsToSearch, HashMap<Integer, Coordinate> stationId2CoordinateMap,
            double[] statValues, HashMap<Integer, Integer> idStationIndexMap ) {

        List<Integer> stationsToUse = new ArrayList<Integer>();
        List<Integer> stationsLeftOver = new ArrayList<Integer>();
        Map<Double, Integer> sortedByDistanceStationsMap = new TreeMap<Double, Integer>();
        for( Integer stId : stationIdsToSearch ) {
            /*
             * check the values of the stations. If there are novalues and there
             * are enough stations left, disable them
             */
            double currentValue = statValues[stationid2StationindexMap.get(stId)];
            if (isNovalue(currentValue)) {
                // if the station has a novalue, jump over it, even without
                // adding the station.
                // out.println("Jump over station: " + stId);
                stationsLeftOver.add(stId);
                continue;
            }

            /*
             * if it gets here, the station has a value. 
             * Put it to the list of stations to be checked 
             * for which is nearer. 
             */
            Coordinate stationCoord = stationId2CoordinateMap.get(stId);
            double distance = basinBaricenterCoordinate.distance(stationCoord);
            sortedByDistanceStationsMap.put(distance, stId);
        }
        Collection<Integer> statIds = sortedByDistanceStationsMap.values();
        Iterator<Integer> iterator = statIds.iterator();
        for( int i = 0; i < statIds.size(); i++ ) {
            if (iterator.hasNext() && i < stationNum) {
                stationsToUse.add(iterator.next());
            } else if (iterator.hasNext() && i >= stationNum) {
                stationsLeftOver.add(iterator.next());
            } else {
                System.out.println("SHOULD THIS EVER HAPPEN???");
                break;
            }
        }

        /*
         * if not enough stations were collected, add also stations that
         * don't have values. Those stations are taken in random way, 
         * since their value won't be considered.  
         */
        if (stationsToUse.size() < stationNum) {
            for( int i = 0; i < stationNum - stationsToUse.size(); i++ ) {
                stationsToUse.add(stationsLeftOver.get(i));
            }
        }

        return stationsToUse;
    }

    /**
     * Fills the elevation and id arrays for the stations, ordering in ascending
     * elevation order.
     * 
     * @throws Exception
     *             in the case the sorting gives problems.
     */
    private void extractFromStationFeatures() throws Exception {
        int stationIdIndex = -1;
        int stationElevIndex = -1;
        pm
                .beginTask(
                        "Filling the elevation and id arrays for the stations, ordering them in ascending elevation order.",
                        stationCoordinates.size());
        for( int i = 0; i < stationCoordinates.size(); i++ ) {
            pm.worked(1);
            SimpleFeature stationF = stationFeatures.get(i);
            Coordinate stationCoord = stationCoordinates.get(i);
            if (stationIdIndex == -1) {
                SimpleFeatureType featureType = stationF.getFeatureType();
                stationIdIndex = featureType.indexOf(stationIdField);
                stationElevIndex = featureType.indexOf(stationElevField);
                if (stationIdIndex == -1) {
                    throw new IllegalArgumentException("Could not find the field: "
                            + stationIdField);
                }
                if (stationElevIndex == -1) {
                    throw new IllegalArgumentException("Could not find the field: "
                            + stationElevField);
                }
            }
            int id = ((Number) stationF.getAttribute(stationIdIndex)).intValue();
            double elev = ((Number) stationF.getAttribute(stationElevIndex)).doubleValue();
            statElev[i] = elev;
            statId[i] = id;
            stationId2CoordinateMap.put(id, stationCoord);
        }
        pm.done();
        // sort
        out.println("Sorting...");
        QuickSortAlgorithm qsA = new QuickSortAlgorithm(pm);
        qsA.sort(statElev, statId);
    }

    private void rangeT( double[] statValues ) {
        // calcola la temperatura massima maxT[j] e minima minT[j] per ogni
        // stazione
        // calcola anche il DTmonth[j], l'escursione termica giornaliera mediata
        // nei 30 giorni precedenti
        // DTday e' variabile ausiliaria

        int stationNum = statValues.length;

        if (cont_min_max == null) {
            cont_min_max = new int[]{0, 0};
            flag_Tmin = new int[stationNum];
            flag_Tmax = new int[stationNum];
            minTempPerStation = new double[stationNum];
            maxTempPerStation = new double[stationNum];
            DTday = new double[stationNum][32];
            DTmonth = new double[stationNum];
        }

        myCalendar.setTime(currentTime);

        int currentHour = myCalendar.get(Calendar.HOUR_OF_DAY);
        int currentMinute = myCalendar.get(Calendar.MINUTE);
        double hh = currentHour + currentMinute / 60.0;

        if (hh_prev != -1) {
            // cicli if da usare nel caso non sia disponibile la temperatura
            // all'istante h_T_min
            // ciclo if per la Tmin a cavallo dell'istante h_T_min
            if (hh >= hourForTmin && hh_prev < hourForTmin) {
                cont_min_max[0] += 1;
                // per ogni stazione di misura di T setto la flag_Tmin a 1
                for( int j = 0; j < stationNum; j++ ) {
                    flag_Tmin[j] = 1;
                }
            }
            // ciclo if per la Tmax a cavallo dell'istante h_T_max
            if (hh >= hourForTmax && hh_prev < hourForTmax) {
                cont_min_max[1] += 1;
                // per ogni stazione di misura di T setto la flag_Tmax a 1
                for( int j = 0; j < stationNum; j++ ) {
                    flag_Tmax[j] = 1;
                }
            }
        }

        // per ogni stazione di misura di T
        for( int j = 0; j < stationNum; j++ ) {

            if (flag_Tmin[j] == 1) { // sono a ridosso della Tmin
                if (!isNovalue(statValues[j])) {
                    minTempPerStation[j] = statValues[j];
                    flag_Tmin[j] = 2;

                    // ADDED
                    if (minTempPerStation[j] >= maxTempPerStation[j])
                        maxTempPerStation[j] = minTempPerStation[j] + DTd_default;

                    // se l'ora di osservazione è fuori dall'intervallo di
                    // tolleranza
                    // setto minT a NODATA
                } else if (hh > hourForTmin + tolleranceHours4Tmin || hh < hourForTmin) {
                    minTempPerStation[j] = JGrassConstants.doubleNovalue;
                    flag_Tmin[j] = 2;
                }
            }

            // come per Tmin lavoro con Tmax
            if (flag_Tmax[j] == 1) {
                if (!isNovalue(statValues[j])) {
                    maxTempPerStation[j] = statValues[j];
                    flag_Tmax[j] = 2;

                    // ADDED
                    if (minTempPerStation[j] >= maxTempPerStation[j])
                        minTempPerStation[j] = maxTempPerStation[j] - DTd_default;

                } else if (hh > hourForTmax + tolleranceHours4Tmax || hh < hourForTmax) {
                    maxTempPerStation[j] = JGrassConstants.doubleNovalue;
                    flag_Tmax[j] = 2;
                }
            }

            // calcolo effettivo del DT giornaliero e aggiornamento della
            // variabile flag_Tmin e max
            if (cont_min_max[0] == cont_min_max[1] && cont_min_max[0] >= 0) {
                if (flag_Tmin[j] == 2 && flag_Tmax[j] == 2) {

                    if (!isNovalue(minTempPerStation[j]) && !isNovalue(maxTempPerStation[j])) {
                        DTday[j][cont_min_max[0]] = maxTempPerStation[j] - minTempPerStation[j];
                    } else {
                        DTday[j][cont_min_max[0]] = DTd_default;
                    }

                    flag_Tmin[j] = 3;
                    flag_Tmax[j] = 3;
                }
            }

        }

        boolean[] hasMinMaxT = new boolean[stationNum];
        boolean atLeastOneStation = false;
        // per ogni stazione di misura di T
        for( int j = 0; j < stationNum; j++ ) {
            if (flag_Tmin[j] == 3) {
                hasMinMaxT[j] = true;
                atLeastOneStation = true;
            }
        }

        if (atLeastOneStation) {
            if (cont_min_max[0] > 30) {
                for( int j = 0; j < stationNum; j++ ) { // j sono le stazioni
                    for( int k = 0; k < 30; k++ ) { // k sono i giorni
                        DTday[j][k] = DTday[j][k + 1];
                    }
                }
                cont_min_max[0] = 30;
                cont_min_max[1] = 30;
            }

            // per ogni stazione di misura di T
            for( int j = 0; j < stationNum; j++ ) {

                if (!hasMinMaxT[j]) {
                    continue;
                }

                // inizializzo a zero la variabile DTmonth
                DTmonth[j] = 0.0;
                int cont = 0;

                for( int k = 0; k <= cont_min_max[0]; k++ ) {
                    if (!isNovalue(DTday[j][k])) {
                        // aggiorno DTmonth con i valori di Dtday e anche il
                        // contatore
                        DTmonth[j] += DTday[j][k];
                        cont += 1;
                    }
                }

                if (cont == 0) {
                    // assegno a NODATA il valore di DTmonth per la stazione
                    // corrente
                    DTmonth[j] = DTm_default;
                } else {
                    // calcolo del Dt medio mensile per la stazione corrente
                    DTmonth[j] /= (double) cont;
                }

                flag_Tmin[j] = 0;
                flag_Tmax[j] = 0;

            }

        } else {
            for( int j = 0; j < stationNum; j++ ) {
                if (DTmonth[j] == 0.0) {
                    DTmonth[j] = DTm_default;
                }
            }
        }

        hh_prev = hh;
    }

    public void addLink( ILink link ) {
        String id = link.getID();

        if (id.equals(valuesID)) {
            valuesLink = link;
        }
        if (id.equals(stationsID)) {
            stationsLink = link;
        }
        if (id.equals(basinsID)) {
            basinsLink = link;
        }
        if (id.equals(outputID)) {
            outputLink = link;
        }
        if (id.equals(altimetryID)) {
            altimetryLink = link;
        }
    }

    public void finish() {
    }

    public IInputExchangeItem getInputExchangeItem( int inputExchangeItemIndex ) {
        if (inputExchangeItemIndex == 0) {
            return valuesInputEI;
        } else if (inputExchangeItemIndex == 1) {
            return stationsInputEI;
        } else if (inputExchangeItemIndex == 2) {
            return basinsInputEI;
        } else if (inputExchangeItemIndex == 3) {
            return altimetryInputEI;
        }
        return null;
    }

    public int getInputExchangeItemCount() {
        return 4;
    }

    public String getModelDescription() {
        return modelParameters;
    }

    public IOutputExchangeItem getOutputExchangeItem( int outputExchangeItemIndex ) {
        return ouputDataOutputEI;
    }

    public int getOutputExchangeItemCount() {
        return 1;
    }

    public void removeLink( String linkID ) {
        if (linkID.equals(valuesID)) {
            valuesLink = null;
        }
        if (linkID.equals(stationsID)) {
            stationsLink = null;
        }
        if (linkID.equals(basinsID)) {
            basinsLink = null;
        }
        if (linkID.equals(outputID)) {
            outputLink = null;
        }
        if (linkID.equals(altimetryID)) {
            altimetryLink = null;
        }
    }

}
